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A nonequilibrium Potts-like model with q absorbing states 
is studied using Monte Carlo simulations. In two dimensions 
and q = 3 the model exhibits a discontinuous transition. For 
the three-dimensional case and q — 2 the model exhibits a 
continuous, transition with /3 = 1 (mean-field). Simulations 
are inconclusive, however, in the two-dimensional case for q = 
2. We suggest that in this case the model is close to or at 
the crossing point of lines separating three different types of 
phase transitions. The proposed phase diagram in the (q, d) 
plane is very similar to that of the equilibrium Potts model. 
In addition, our simulations confirm field-theory prediction 
that in two dimensions a branching-annihilating random walk 
model without parity conservation belongs to the directed 
percolation universality class. 



I. INTRODUCTION 

One of the main achievements of statistical physics 
during the past decades is the understanding of the uni- 
versal properties of systems near an equilibrium second 
order phase transition. Universality classes, character- 
ized by a small number of parameters allow us to under- 
stand why different systems have the same critical prop- 
erties 0. 

One important model in equilibrium statistical physics 
is the so-called q-states Potts model g], which can be 
used to describe a large class of physical systems. Its rich 
critical behaviour is to a large extent understood. It is 
known that the important parameters which determine 
its critical behaviour are its dimensionality d and the 
degeneracy of its ground state q. 

The relatively complete understanding of equilibrium 
phase transitions has yet no nonequilibrium counterpart. 
However, it is becoming evident that some analogies be- 
tween equilibrium and nonequilibrium systems could be 
made. 

For example, a classification into universality classes is 
particularly evident for one-dimensional nonequilibrium 
models with absorbing states ||, which exhibit a phase 
transition between active and absorbing phase in their 
stationary state. A prime example of a universality class 
is the directed percolation one, which is typical to models 
with a single absorbing state [Q. The so-called parity- 
conserving (PC) universality class is typical for models 



with certain conservation laws or symmetries . 

However, it is necessary to investigate nonequilibrium 
models in higher dimension to have a complete charac- 
terization of the possible universal properties. 

Analytical approach to high dimensional problems is 
based mainly on field-theory methods. Recently, inter- 
esting results were obtained along this line by Cardy and 
Tiiuber || who clarified the role of parity conservation 
and also discovered some new universality classes. How- 
ever, this technique is applicable only to certain parti- 
cle systems |[o| and a large class of two- and higher- 
dimensional models with absorbing states cannot be 
treated within such a method. Accordingly, our under- 
standing of higher dimensional nonequilibrium models is 
rather limited. As for continuous phase transitions some 
models with single absorbing state were shown to fall 
into the DP universality class Jll[]. Results were also re- 
ported for a two-dimensional model with two absorbing 
states H . These models are often motivated by problems 
of surface catalysis In some models with absorbing 
states discontinuous phase transitions are also known to 
occur ©10. 

A subclass of models with absorbing states are the 
so-called branching annihilating random walk (BARW) 
ones. In BARW models each particle can react (annihi- 
late, branch, diffuse, etc.) according to prescribed rules. 
It turns out that BARW models for d > 1 are much differ- 
ent from the above mentioned (surface-catalysis) models. 
Some of the results obtained from field-theory method |]] 
have been confirmed using numerical methods. For ex- 
ample Szabo and Santos confirmed the existence of loga- 
rithmic corrections in two-dimensional parity-conserving 
BARW models pq| . However, for parity non-conserving 
particle systems Monte Carlo simulations J B[ seem to be 
in disagreement with field-theory results |9 1 . 

A characterization of the rich critical behaviour en- 
countered in models with absorbing states is clearly an 
important issue. However, lacking a sound theoretical 
basis, it is by no means obvious which parameters are 
relevant for such a classification. On general grounds 
one expects that the dimensionality d of the system is a 
relevant parameter. Moreover, based on the information 
coming from the one- dimensional case and from equi- 
librium one expects that the number of absorbing states 
q is another relevant parameter. Note however that, in 
some cases, details of the dynamics might also change 
the critical properties, even without affecting the sym- 
metricity of absorbing states p6). Nevertheless, we can 
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expect that these cases are accidental rather than generic 
(as it is the case in equilibrium phase transitions when a 
marginal scaling field is present ]l7p ). Keeping in mind 
the above objections, and lacking a better candidate, we 
consider q as a classification parameter |Q . 

In the present paper we study a recently introduced 
noncquilibrium Potts model with q absorbing states. Our 
numerical results for d > 2 can be summarized as follows: 
(i) For d = 2 and q = 3 the models has a discontinuous 
transition, (ii) For d = 3 and q = 2 the model exhibits 
a continuous, transition with an order parameter critical 
exponent f3 — I (mean-field), (iii) Simulations are incon- 
clusive, however, in the d = 2 and q = 2 case, and we 
suggest that in this case the model is close to or at the 
crossing point of lines separating three different types 
of phase transitions. Our results together with the al- 
ready accumulated knowledge, prompted us to partition 
the (q, d) plane into three regions of different phase tran- 
sitions: mean-field, non-mean-field, and discontinuous. 
Somewhat surprisingly, the topology of such a partition 
is the same as for the equilibrium Potts model. Although 
nonequilibrium systems, and in particular models with 
absorbing states, are usually regarded as very much dif- 
ferent from equilibrium systems our work shows, how- 
ever, that despite some differences there are also some 
qualitative similarities. 

In addition, we performed simulations of a d = 2 
BARW model without parity conservation. Our results 
confirm the field-theory predictions according to which 
the critical behavior belongs to to the directed percola- 
tion universality class. 

II. NONEQUILIBRIUM POTTS MODEL 

Before presenting our model let us recall basic proper- 
ties of the equilibrium Potts model. First we assign with 
a lattice site i a g-state variable oi — 0, 1, q— 1. Next, 
we define the energy of this model through the Hamilto- 
nian: 

(ij) 

where summation is over pairs which are usually 

nearest neighbours and 5 is the Kronecker delta function. 
This equilibrium statistical mechanics model was studied 
using many different analytical and numerical methods 
and is a rich source of the information about phase tran- 
sitions and critical phenomena Q . 

To study the equilibrium Potts model using Monte 
Carlo simulations one constructs a stochastic Markov 
process with suitably chosen transition rates. One of the 
possible choices corresponds to the so-called Metropolis 
algorithm. In this algorithm one looks at the energy dif- 
ference AE between the final and initial configuration 



and accept the move with probability min 
where T is temperature. 

A nonequilibrium Potts model having q adsorbing 
states can be obtained by making the following transfor- 
mation in the Metropolis dynamics p6| : when all neigh- 
bours of a given site are in the same state as this site, 
then this site cannot change its state (at least until one 
of its neighbours is changed). Let us notice that any 
of q ground states of the equilibrium Potts model (jl]) 
is an absorbing state of our noncquilibrium Potts model. 
The one-dimensional version of this nonequilibrium Potts 
model has already been examined |l6| . In addition to re- 
covering the expected critical behaviour for q — 2 and 3 it 
was found that certain additional modification of its dy- 
namics might affect the critical behaviour. We shall not 
be concerned with such a variant in the present paper. 

Since the dynamics of our models is obtained from a 
modification of the Metropolis algorithm of an equilib- 
rium system, transition probabilities are parametrized by 
temperature-like quantity T. Strictly speaking, for our 
model the ordinary temperature cannot be defined. Nev- 
ertheless, we will refer to this quantity as temperature. 

Of course in the realm of nonequilibrium systems there 
are also other models than those with absorbing states. 
One of the important questions is whether under certain 
conditions nonequilibrium systems might be mapped, at 
least at the coarse-grained level, into equilibrium ones. 
Some aspects of this problem were studied by Grinstein et 
al. . In a class of models studied by them it is impor- 
tant that all transition probabilities are strictly greater 
than zero. In models with absorbing states this require- 
ment is clearly violated. 

III. MONTE CARLO SIMULATIONS AND 
RESULTS 

To study the properties of our model we made exten- 
sive Monte-Carlo simulations. A natural characteristic of 
models with absorbing states is the steady-state density 
of active sites p. A given site i is active when at least 
one of its neighbours is in a state different than i. Oth- 
erwise the site i is called nonactive. Upon approaching a 
critical point p develops a power-law singularity charac- 
terized by the exponent /?. At the first-order transition 
p has a jump. In addition, we also looked at its time 
dependence p(t). In the active phase p(t) converges to 
the positive value, while at criticality p(t) usually has a 
power-law decay p ~ t~ s . In the absorbing phase p very 
often decays faster than the power law, however, in some 
cases (also those studied in the present paper) a power 
law behaviour is seen, but with a different exponent than 
for the critical decay. 

Moreover, we used the so-called dynamic Monte Carlo 
method where one sets the system in the absorbing state 
with activity only locally initiated and monitor some 
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stochastic properties of runs [M ■ O ne of the most fre- 
quently used characteristics is the survival probability 
P(t) that activity survives at least until time t and the 
average number of active sites N(t) (to calculate N(i) 
we average over all runs). At criticality P{t) and N(t) 
are expected to have power-law decay: P(t) ~ t~ 5 and 
N(t) ~ t v . (For some models 6 = 5', but exceptions from 
this relation are also known [0]). Our simulations were 
made for various system sizes and we ensured that the 
system was large enough so that presented below results 
are size independent. 



d = 2, q = 3 

Simulations for d = 2 models were performed on the 
square lattice. The temperature dependence of the den- 
sity of active sites p is shown in Fig. |]. For T < 1.237 it is 
virtually impossible to reach an active steady state value 
of p, which suggests that the model undergoes a discon- 
tinuous phase transitions. Such a scenario is confirmed 
by the data in Figs. ||{|. In Fig. ||one can see that upon 
approaching the transition point p(t) develops a longer 
and longer plateau. At the transition point, which we 
locate around T = 1.237 the density p(t) has basically 
an infinitely-long plateau. Such a behaviour is a clear 
indication of the discontinuous nature of the transition. 
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FIG. 1. The density of particles p as a function of temper- 
ature T for the d — 2 nonequilibrium Potts model. Close to 
the transitions we used the lattice of the linear size L = 600. 

Let us notice that below the transition temperature 
p(t) decays as t^ 1 ^ 2 . A simple scaling argument can be 
used to show that such a behaviour is related with the 
average domain size growth I ~ t 1 / 2 . Indeed, let us con- 
sider the d = 2 system of the linear size L. It contains 
(j) 2 subdomains of the linear size I and thus the total 
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density scales as jj-j = l~ L . Assuming now that I in- 
creases as t 1 / 2 we obtain that p ~ t^ 1 / 2 . The character- 
istic length which increases as t 1 / 2 is typical for coarsen- 
ing in the broken-symmetry phase of the original Potts 
model (2|j2^]. It also appears in the evolution of one- 
dimensional nonequilibrium systems which belong to the 
parity-conserving universality class. 
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FIG. 2. The time dependence of the density p(t) for the 
d = 2 model and (from top) T =1.3, 1.25, 1.237, 1.2365, 
1.236, 1.235, 1.23, 1.2, and 1.1 (L = 500). Each line is an 
average of 100 independent runs which starts from random 
initial configurations. The dotted straight line has a slope 
corresponding to 8 = 0.5. 
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FIG. 3. The time dependence of the survival probability 
P(t) for the d = 2 model and (from top) T =1.25, 1.245, 
1.242, 1.241, 1.24, 1.239, 1.238, 1.237, 1.2365, 1.236, and 1.23 
(L = 500). Lines for T = 1.2365, 1.237, and 1.238 are ob- 
tained averaging over 10 8 independent runs. 
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FIG. 4. The time dependence of the number of active sites 
N(t) for the d = 2, g = 3 model and (from top) T =1.25, 
1.245, 1.242, 1.241, 1.24, 1.239, 1.238, 1.237, 1.2365, 1.236, 
and 1.23 (L = 500). Lines for T = 1.2365, 1.237, and 1.238 
are obtained averaging over 10 8 independent runs. 

Results of the dynamic Monte Carlo also support the 
discontinuous nature of the transition. Indeed, both for 
P(t) (Fig. |) and N(t) (Fig. |) the data are systemat- 
ically bending and no clear power-law behaviour is ob- 
served. 

Thus, as a summary, for the d = 2, q = 3 case our 
model undergoes a discontinuous phase transition. 



d = 2, q = 2 



where m is a certain constant. Application of this method 
to our data and m = 5 is shown in Fig. ||. At criticality 
the data converge (with some scattering) to the value 
0.06 but such a small value means that the possibility 
5 = still cannot be excluded. 

In our opinion, from the steady-state and time- 
dependent measurements of p the most likely possibility 
is that the transition is of first order and the plateau will 
develop (in Fig. 0) but only at a larger time scale. 




FIG. 5. The time dependence of the density p(t) for the 
d = 2 model and (from top) T =1.77, 1.76, 1.7585, 1.758, 
1.757, 1.75, 1.74, and 1.71 (L = 500). Each line is an aver- 
age of 100 independent runs which starts from random initial 
configurations. The dotted straight line has a slope corre- 
sponding to 8 = 0.5. 



In the q — 2 case the density of active sites is a similar 
function of temperature as in the q = 3 case (Fig. [l]). 
Let us notice however, that now the jump in p is almost 
twice smaller than previously. We should be also aware 
of the fact that an observed jump might be a finite-size 
effect. (Large fluctuations might drive the system into 
an absorbing state even though temperature is above the 
critical temperature). As further results show, it is rather 
difficult to clarify the nature of the transition in this case. 
First, let us notice that time dependence of pit") does not 
develop a clear plateau as in the q = 3 case. (We estimate 
the critical point in this case as T c = 1.7585(3)). But 
neither there is a pronounced power-law behaviour seen 
in Fig. [|. If pit) does decay as t~ s at criticality then 
5' is very small (<5' = 0.07(2)), which suggests that the 
true exponent 5' might be equal to zero. Sometimes to 
improve the estimation of 6 one can use the so-called local 
slopes method ||. In this method 5 is calculated from 
the equation 
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FIG. 6. The exponent 5(t) as a function of 1/t for the 
d = q = 2 model and (from top) T =1.77, 1.7585, 1.757, 1.75, 
1.74 and 1.71 (L = 500). Each line is an average of about 100 
independent runs. 

However, the dynamic Monte Carlo results shed some 
doubts on such an interpretation. Indeed, in Fig. ^ we 
can see that P(t) has a clear power-law decay for at least 
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three decades in time with the exponent 5' — 0.90(2). In 
addition, N(t) seems to remain constant at the transi- 
tion point which suggest that -q — 0. To support our dy- 
namic Monte Carlo results let us notice that Hinrichsen 
reported basically the same values for dynamical expo- 
nents for another two-dimensional model with q — 2 m. 



some dynamical power-law characteristics, although ex- 
otic at first sight, cannot be ruled out. We will return to 
this problem in the next section. 

Let us also notice (Fig. ||) that, similarly to the q = 3 
case, for T < T c the density p(t) seems to decay in time 
asi- 1 ^. 
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FIG. 7. The time dependence of the survival probability 
P(t) for the d = 2 model and (from top) T =1.77, 1.76, 1.759, 
1.7585, 1.758, 1.757, 1.755, 1.74, and 1.73 (L = 500). Each 
line is an average of about 10 6 independent runs. 
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FIG. 8. The time dependence of the number of active sites 
N(t) for the d = 2, q = 2 model and (from top) T =1.77, 1.76, 
1.759, 1.7585, 1.758, 1.757, 1.755, 1.74, and 1.73 (L = 500). 
Each line is an average of about 10 6 independent runs. 

The fact 5' is close to unity and rj = prompted him 
to suggest that the model exhibit a mean field behaviour. 
In such a case, however, one should have 5 = 1, which is 
clearly in contradiction with the behaviour of pit) which 
we observe. A power-law behaviour of P(t) and N(t) is 
typical for continuous transitions which seems to be in 
conflict with the behaviour of p. However, the possibility 
that we have a discontinuous transition accompanied by 



d = 3, q = 2 

As a last case in this section, we consider our model on 
the simple cubic lattice and for q = 2. Of course, simula- 
tions for three-dimensional systems are very demanding. 
Therefore we were not able to perform detailed time- 
dependent simulations nor the dynamic Monte Carlo. 
However, the steady-state density p as measured for 
L = 40 and 60 is basically size independent (Fig. ||). 
Nearly linear behaviour of p in the vicinity of the tran- 
sition suggests that in this case (3 — 1, which indicate a 
mean-field nature of the transition. 
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FIG. 9. The density of particles p as a function of temper- 
ature T for the three-dimensional q — 2 nonequilibrium Potts 
model with L = 60(D) and 40(+). 

(q,d) phase diagram 

In this section we sketch the overall behaviour of our 
model in the (q, d) plane. Our proposal (Fig. [To]) is based 
on the already accumulated knowledge, presented above 
results, and a minimalistic assumption that the resulting 
phase diagram should not be too complicated. Essen- 
tially, we suggest that the (q, d) plane can be divided 
into three regions of different types of phase transitions: 
(i) non-mean-field, (ii) mean-field, and (iii) first-order. 
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FIG. 10. The location of various types of phase transitions 
for d— dimensional models with q absorbing states. Models for 
which numerical results are reported in this paper are denoted 
as stars. 

Some comments are in order. As it is defined only- 
through dynamical rules, our model is trivial for q = 1. 
However, the q = 1 case corresponds to models with a 
single absorbing state (directed percolation, contact pro- 
cess) It is known that the critical dimension d c in this 
case equals 4 and for d smaller/larger than d c we have a 
continuous non-mean- field/mean- field transition. Firmly 
established are also results along the d = 1 line. It is 
known that q = 2 case corresponds typically to the PC 
universality class. For q > 3 one expects |||23| that the 
model typically belongs to the same universality class 
as a A-BARW model studied by Tauber and Cardy ||, 
although, under more restrictive dynamics the PC crit- 
icality might also appear |16[) . In any case, the criti- 
cal behaviour is non-mean- field. Taking into account the 
mean-field behavior obtained in the case d = 3, q = 2, 
we assume that this case falls into the same region as DP 
above critical dimension. On the other hand, the dis- 
continuous transition in the d = 2, q = 3 case implies 
the existence of the third region. Actually we performed 
some simulations also for d = 2, q = 4 case. Basically 
the behaviour is similar to the q — 3 case except that the 
jump is larger and the first-order character of the tran- 
sition is even more transparent. We expect that such a 
behaviour persists for all q > 3. Provided that for d > 2 
there are no qualitative changes, the diagram must have 
a structure as shown in Fig. 

It is thus clear that our most difficult case d = q = 2 
is located somewhere close to the point where all three 
regions meet. We have no strong arguments to locate this 
point exactly at the crossing point, but such a location 
would certainly explain unusual behaviour as seen in our 
Monte Carlo simulations. We do not exclude, however, 
the possibility that d = q = 2 case is off the crossing point 
but somewhere close to it, which would still explain the 



numerical difficulties in this case. 

The resulting diagram looks very similar to the dia- 
gram of the equilibrium Potts model . In the equilib- 
rium model the non-mean-field part is basically shifted 
by two upwards (critical dimensions for ordinary perco- 
lation and the Ising model are 6 and 4, respectively). 
However, the diagram for the equilibrium Potts model 
is much more meaningful than in our case. Indeed, in 
the Fortuin-Kasteleyn representation, the Potts model is 
well defined for any, even non-integer q, which justifies 
continuous lines on its phase diagram. In our case, we do 
not have a representation of our model with continuous 
q. The key property that would be required for such a 
representation is the existence of an analog of the parti- 
tion function. Such a quantity exists for nonequilibrium 
systems only in very special cases Provided that a 
certain partition function exists for our model (no mat- 
ter how complicated), and that using this function one 
can find a corresponding Fortuin-Kasteleyn representa- 
tion (no matter how complicated), continuous lines in 
our diagram would be meaningful. 

Of course, the presented diagram does not encompass 
all models with absorbing states. It is well known that for 
example there are d = 2, q = 1 models with first-order 
transitions . But as we already mentioned, similar situ- 
ation occurs for equilibrium systems where certain factors 
(anisotropies, additional interactions etc.) might change 
more generic behaviour. The presented diagram is valid 
only for the presented Potts model and its applicability 
to other systems requires additional examination. 



IV. PARITY-NON CONSERVING BARW MODEL 
IN TWO DIMENSIONS 

In the present section we examine the two-dimensional 
BARW model without parity conservation. In this model 
particles are located on sites of a square lattice. In addi- 
tion to diffusion, which takes place at the rate p, particles 
can branch, at the rate 1 — p, according to the following 
reaction: 



X -> 2X, 



(3) 



where the offspring particle is placed on the randomly 
chosen nearest neighbour of a parent particle. Moreover, 
two particles which happen to be placed at the same site 
annihilate instantaneously 



2X -> 



(4) 



This model was already examined by Takayasu and 
Tretyakov || . They suggested that the model undergoes 
a continuous transition around p = 0.85 and the den- 
sity of particles decays linearly at the transition point 
(P = 1). This result contradicts more recent field-theory 
approach which suggests that in this case the model 
should belong to the DP universality class M. 
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FIG. 11. The density of particles p as a function of the 
diffusion rate p for the two-dimensional BARW model with 
L = 1000(D) and 2000(+). 

Since such a disagreement requires an explanation, 
we performed Monte Carlo simulations of this model. 
Our system size was much larger than in Takayasu and 
Tretyakov simulations and we approached much closer to 
the critical point. One can see (Fig. [ll]) that although 
around p = 0.8 (which was the largest value of p simu- 
lated by Takayasu and Tretyakov) the density seems to 
decay linearly it has a pronounced bending close to the 
transition point. To obtain more accurate estimation of 
the critical point wc examined the time dependence of 
p(t) (see Fig. |lj). From these analysis we obtain the 
following estimation of the critical point p c = 0.8237(5). 
One can also see that at criticality p(t) has a power-law 
decay with the exponent close to the DP value S = 0.451. 
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Having the critical point we can estimate exponent (3 
and the corresponding data are shown in Fig. [l^. The 
least-square fit gives f3 = 0.60(3) which is certainly com- 
patible with the DP value 0.584(4). 
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FIG. 13. The density of particles p(t) as a function of 
p c — p for the two-dimensional BARW model (p c — 0.8237). 
Only results for L = 2000 runs are shown here. The linear 
fit (dotted line) obtained using the least-square method has 
a slope corresponding to /3 = 0.6. 

To summarize this section, our results confirm the 
field-theory prediction that BARW models without par- 
ity conservation belong to the DP universality class. 



V. CONCLUSIONS 

In the present paper we examined d-dimensional 
noncquilibrium models with q absorbing states. As our 
main result we obtained the diagram shown in Fig. [l^. 
Interestingly, this diagram bears some similarity to the 
diagram of equilibrium Potts model. In addition, we clar- 
ified the nature of the phase transition in the d = 2 
BARW model without parity conservation. Together 
with the work of Szabo and Santos Jl5| for the parity 
conserving case, it confirms predictions of the field-theory 
for d = 2 BARW models by Cardy and Tauber §. 

Although it requires considerable numerical efforts, it 
would be desirable to clarify the behaviour of the d = q = 
2 model. Our results are inconclusive in this case, but a 
possibility that an interesting critical behaviour could be 
found should motivate further study. 



FIG. 12. The time dependence of density of particles 
p(t) for the two-dimensional BARW model. Simulations were 
made for (from top) p =0.8, 0.81, 0.82, 0.822, 0.823, 0.8237, 
0.824, and 0.825. The dotted line has a slope corresponding 
to the DP value S = 0.451. 
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